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solution of one-dimensional Euler 
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term 
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Abstract 


In this paper a novel Godunov-type finite volume technique is presented 
for the solution of one-dimensional Euler equations. The numerical scheme 
defined herein in is well-balanced and approximates the solution by propa- 
gating a set of jump discontinuities from each Riemann cell interface. The 
corresponding source terms are then treated within the flux-differencing of 
the finite volume computational cells. First, the capability of the numerical 
solver under gravitational source term is examined and the results are val- 
idated with reference solution and higher-order WENO scheme. Then, the 
well-balanced property of the scheme for the steady-state is tested and finally 
the proposed method is employed for the modeling small and large amplitude 
perturbation imposed to the polytropic atmosphere. It is found out that the 
defined well-balanced solver provides sensible prediction for all of the given 
test cases. 


Keywords: Wave propagation algorithm, Flux wave formula, Riemann 
solver, Well-Balanced, Euler equations. 


1. Introduction 


The Euler equations with gravitational source terms have been extensively 
used in many scientific aspects such as aerospace, astrophysics, and shock 
tube problems. Prediction models should be able to capture sharp gradi- 
ent shocks as well as rarefaction waves that appear within the solution in 
particular with existence of source terms. Generally, two different class of 
finite volume methods have been used to model Euler equations in recent 
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years mainly reviewed by LeVeque [6),{@jand [16]. The first method is upwind 
method, which basically uses the Godunov scheme. Despite the complexity 
of upwind schemes in particular in dealing with associated Jacobian matrix 
they provide very accurate results for the shock capturing problems [I]. 
The second approach is central schemes, which applies Lax—Friedrichs or 
Lax-Wendroff methods. However, they produce rather high diffusion, which 
eventually affects the methods stability unless a refined mesh is used [Z]. 


Previously, significant attentions have been paid for the solution of the 
Euler equations with the gravitational source terms mostly based on the finite 
volume methods. LeVeque and Bale [&] have developed a quasi-state wave 
propagation algorithm for the solution of the Euler equations. In another 
work Botta et al. [2] defined a well-balanced finite-volume methods, which 
maintain certain class of steady states for the nearly hydrostatic flows. A well- 
balanced approach on the basis of the gas-kinetic scheme has been proposed 
by Luo et al. [10] for an isolated gravitational hydrodynamic system. Kappeli 
and Mishra |5] have introduced a second-order well-balanced finite volume 
scheme for the isentropic hydrostatic equilibrium. Chandrashekar and Klin- 
genberg [8] implemented a well-balanced second-order Godunov-type finite 
volume method for the Euler equations with gravitation. More recently, Li 
and Xing designed a high-order well-balanced finite volume WENO (weighted 
essentially oscillatory) scheme for the Euler equations with gravitational field, 
which preserves both isothermal equilibrium and the polytropic hydrostatic 
balance state. 


The main purpose of this work is to develop a version of Godunov-type 
wave propagation algorithm for the solution of one dimensional Euler equa- 
tions. The proposed method extends the a modified flux-wave solution pro- 
vided in [11,2] for the shallow water equations (SWEs) to the one dimen- 
sional Euler equations. This approach is well-balanced and treats any source 
terms within the flux-differencing of the finite-volume computational cells. 
Additionally, the defined numerical scheme utilizes the advantage of com- 
bination both approximate and exact Riemann speeds, which enables the 
method to avoid non-negative pressure fields for the Euler equations. To 
the best of author’s knowledge no development of the modified flux-wave 
approach defined in [11,2] is used for the solution of the one dimensional 
Euler equations with gravitational source term. The rest of this paper is 
organized as follows: In the next section the mathematical equations for the 
one dimensional Euler equations consisting of gravitational source terms are 
provided. Then in the second section, the wave propagation algorithm with 
both first-order and high-resolution accurate terms is expressed. In the third 
section, the flux-wave formula comprising different choice of wave speeds for 
the one dimensional Euler equations is described. Fourth section states the 
validation of the introduced numerical method with the reference solutions 
and other numerical results available in literature. Finally, the paper ends 
with the summary of numerical results and conclusions of findings. 
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2. Governing equations 


The one dimensional Euler equations including source terms can take the 
conservation law form as 


U;t+F(U), =S, 
p pu 0 
U= pu , F(U) = pu? +P ,S= —pox ? 
E (E+ P)u —puds 


where U is the vector of unknowns, F(U) is the flux-term, S shows cor- 
responding source term, p is density,u denotes the particle velocity, P is 
pressure,@ is an time independent gravitational potential, and finally EF is 
the total energy, which can be obtained as 


PA 
B= ——+ =p? 
a 


where y represents the ratio of specific heat. The relevant eigenvalues and 
eigenvectors for the defined system of equations are expressed as 


Ay =u-—c, A2 =u, A3 =Utc. 


1 1 1 
ry = | u-c|,r2= U ,7T3= | utc, 
¢—uc u? /2 ¢ + uc 


where in the above equations c and ¢ are called sound speed and the total 
specific enthalphy, respectively, which can be computed as 


yP E+P 
c= = : 
p p 


In order to solve the above system of equation the wave propagation algorithm 
described in the next section is used. 


3. Wave propagation algorithm 


The one dimensional Godunov-type wave propagation algorithm can be given 


as (6,17) 


At 


x 


= At =n ~n 
(ATAU ;_1/2 + ATAU 4414/2) — 5 (Fisspe = F472), 


Ul — U" = 
Ax 
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where U"*" is the vector of unknowns at the next time step, ATAU ;_1 /2 and 


A’ AU;_1/2 provide the right- and left-going fluctuations, and finally F*., 
shows the second-order correction terms required to obtain high-resolution 
scheme, which can be used with different choice of limiters. If F = 0, the 
first-order Godunov-type wave propagation algorithm is achieved. The right 
and left-going fluctuations are then computed, using the following formula- 
tions 


A*AU{_1/2 = > See 1pos AAU ;_1/2 = S- 41/29 


k:Aj_1/2>0 ki:Aj_1s2<0 


where €;,;~1/2 is called the kth flux-wave propagating from cell interface 
i — 1/2 and can be obtained by multiplying particular coefficients into its 
corresponding eigenvector, say, 


€4 1-1/2 = Pri—-1/2V k,i—1/2- 


4. Flux-wave formula 


The flux-wave approach has been first introduced by [[]] for the acoustic 
problem. This approach has been later modified by Mahdizadeh et al. [1,12] 
for the SWEs with modeling wet/dry front abilities. In this section this 
modified version of flux-wave approach is extended for the solution of one 
dimensional Euler equations with gravitational source term .The original flux- 
wave formula including the treatment of source term can take the form [I] 


My 
F(U;) — F(Ui-1) — Si-1p2Aa = ba €ki—1/2> 

k=1 
where F(U;) and F(U;_1) are the fluxes at the left and right side of cell in- 
terface i—1/2 and M,, denotes the number of waves, which for the prescribed 
Euler equations is equal to three and Az implies finite-volume cell length. 
To expand the flux-wave approach for the one dimensional Euler equations 
it is only required that the differences between neighboring fluxes minus the 
source terms is equalized with the summation of the relevant fluxes . This 
can be accomplished as 

Woodward-—Coella 


pity Pi-1Ui-1 0 
pitti +P, || pi-rt?@y+Pi-1 |-Ax | —pi(giti — 6i)/Ax | (1) 
(EB; + Pu; (Bi a) —piti(dis1 — bi)/Ax 
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1 1 1 
=f Uy =c + Bo Uy + B3 Ui + G : 
Ci — Wis ii; /2 Ci + Hcy 


where u and é are the velocity and total specific enthalphy again, which can 
be obtained through the combination of exact and approximate Riemann 
wave speeds fully explained in [D1], where the approximate Riemann solver 
utilized herein is based upon the Roe solver [14] . For the Euler equations 
the approximate velocity and total specific enthalphy can be givens as 


VPi-1wi-1t Jpiti x VS Pi-1G-1 + JPG 
GPa tle JP tale. * 


and the sound speed ¢; can take the form 


“= 


a = 7-1 - 1/288). 


The systems mentioned in equation () can be rewritten as 


i a 4 3, o 
w-G Uw WtE Bo} = | do], (2) 
G, — tiie, 03/2 G + tic: | | Bs 63 


where 61, 69, and 63 are 


1 pit; — Pi—1Ui-1 
62} = | (pitt; + Pi) — (pi-10?_, + Pi-1) + pi(dia1 — 4%) 
63 (Be Pits = Yai > Pea Git + pita Gir = Oe) 


By solving the linear system provided in equation (Q), the coefficients 61, 32, 
and {$3 are computed, which can be later used to calculate the flux-wave 
€ 1-1/2 Tequired for obtaining the left- and right-going fluctuations for the 
first-order Godunov-type wave propagation algorithm. It should be stressed 
for the solution of linear system given above the LU decomposition algorithm 
with partial pivoting defined in [3] at each time-step. 


4.1 Stability conditions 


To ensure the method’s stability, the Courant—Friedrichs—Lewy condition 
(CFL) [4] similar to the SWEs is used. This condition can be given as the 
following equation for the one dimensional Euler equations: 


CFL = 
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where A = maxz(Aj, 2,3) is the maximum amount of wave speeds, where 
each wave speed is obtained as 


M = tj — &, Ao = tH, Ag = tH + &. 
Generally, the wave propagation algorithm uses values of CFL number close 


to one, which ultimately reduces the computational setup time compared to 
other Riemann solvers provided in literature. 


5. Numerical results 


In order to examine the effectiveness of modified flux-wave approach (MF W) 
for the solution of the one dimensional Euler equations, in this section sev- 
eral numerical test cases are provided. The suitability of the proposed MFW 
approach in dealing with the gravitational source terms is first investigated, 
and the obtained numerical results are compared with both reference solu- 
tions and available numerical data borrowed from literature. Then, the well- 
balanced property of the defined method with the existence gravitational 
field for the isothermal equilibrium is verified, and the calculated numerical 
results are validated with the nonwell-balanced WENO scheme. Finally, the 
capability of the proposed well-balanced model in capturing small and large 
amount of amplitude perturbation is tested. It should be expressed that the 
MFW model defined herein was solved, using an in-house FORTRAN code 
on an Intel Core (i7-4790) 3.6 GHz processor with 16GB of RAM. Addition- 
ally, the solver employs high-resolution wave propagation algorithm based on 
the monotonized centered (MC) limiter. The number of computational finite 
volume cells is defined separately for each test-case. 


5.1 Shock tube problem with gravitational field 


The purpose of this test case is to assess the performance of the proposed 
numerical solver for the solution of the one dimensional Euler equations with 
the existence of gravitational source term. For this problem the Sod test case 
is considered with the initial condition as 


(pa PR) = (1, 0, 1) if x < Os 

(0.125,0,1) otherwise, 
where p°,u°, and P® are again the initial data for the defined Riemann prob- 
lem at time t = 0. As the boundary conditions the extrapolation boundary 
conditions were imposed for both left and right boundaries. In order to con- 
sider the effect of gravitational source term a constant exhibit gravitational 
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field 6, = 1 is used within the source term. The computation is implemented 
until time t = 0.2s with 100 uniform cells and the CFL number equal to 0.9. 
Figure 2 shows numerical predictions for the density, pressure, energy and ve- 
locity obtained based on the MFW approach. The validations are performed 
with the reference solution achieved using 2000 numerical cell and also with 
a novel higher-order WENO scheme defined by [9] with the same computa- 
tional volumes. For the density field a left-going rarefaction wave along with 
the contact discontinuity and shock wave are created within solution and the 
obtained results are in qualitative agreement with both reference solution and 
the WENO approach. Additionally, due to the effect of gravitational source 
term the density is pulled upward for the left boundary. 


Figure [I(b) demonstrates the numerical results for the pressure field. As 
can be observed, a left-going rarefaction waves as well as right-moving shock 
are appeared within solution and again the agreement between the obtained 
numerical results and the reference solution are quite well. Figures [(c) and 
M(d) exhibit the numerical results for the energy and velocity, respectively. 
In terms of the energy field a rather similar shape to the pressure is seen. For 
the velocity distributions a negative velocities in some regions are appeared, 
which is mainly caused by gravitational term. 


5.2 One dimensional gas-falling into fixed external 
potential 


This test case was originally introduced in [15,[18] and utilized to verify the 
well-balanced property of the proposed MFW approach for the isothermal 
equilibrium within gravitational field. The initial conditions of the problem 
can be given as 


ir) u=0, and p = RT (-¥r), (3) 


where gravitational field is expressed as 


P = Poexp (- 
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Figure 1: Numerical solution of the shock tube problem with gravitational source term 
performed at t = 0.2s. (a) Density, (b) pressure,(c) energy,(d) velocity 


whereL is the computational domain set equal to 64 and other parame- 
ters are: pp = 1, y = 5/3, T = 0.6866 and ¢o = .02. The solution was then 
performed in double precision using 200 and 400 computational cells until 
time t=50s, where the steady-state condition is reached. Table [1] demon- 
strates the Euclidean norm achieved for the conserved variables p, pu and E 
at the steady-state condition. As can be seen for all variables relatively small 
amount of error is obtained, which clearly state that the initial conditions 
have been preserved during the simulations. 

In order to compare the performance of the defined MFW approach with 
the unbalanced scheme a small perturbation was imposed into the steady- 
state solution provided in (B).Therefore the initial condition related to pres- 
sure becomes 


= AE (-z) + 0.001lexp (—10(x — 32)*). 
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Table 1: Euclidean error norm computed based on the MFW approach for the steady- 
state 


Number of Cells p pu E 
200 2.6493E-04 9.7563E-05 4.3708E-04 
400 6.8740E-05 2.2979E-05 8.9109E-05 


The simulation was implemented up to time t = 179980.53s, which is much 
larger than the reference sound crossing time approximated as T = 120s. Fig- 
ure illustrates the numerical results obtained based upon the second-order 
MFW approach in comparison with non-well-balanced WENO scheme bor- 
rowed from [9]. As can be observed the initial condition for the velocity pro- 
file was maintained during computation in contrast to the non-well-balanced 
scheme, which was not able to properly neutralize the effect of gravitational 
source terms with the flux gradient for the isothermal equilibrium state. 


5.3 Small and large amplitude wave propagation 


The final test case employed herein was introduced in [5] and investigates the 
performance of the well-balanced MFW approach in dealing with small and 
large disturbances for the polytropic atmosphere in the gravitational field 
(x) = gx. In doing so, the polytropic steady state solution is defined as 


1y-1 


ie (oo : mie) Oe tia Venta = eptay™. 


with the constants: g = 1, y = 5/3, po = 1, po = 1, and Ko = po/(pQ). 
The computational domain of the polytropic atmosphere was set equal to 
[0,2] and the following periodic velocity perturbation was imposed at the 
bottom boundary 


u(0,t) = Asin(4zt), (4) 
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Figure 2: The numerical results achieved with the MFW method along with the 
non well-balanced WENO scheme for the perturbation problem given in section 5.2 at 
t = 179.980.53s. (a) Density, (b) pressure,(c) velocity 


In the first numerical study the amount of A in above equation was chosen 
to 10~®which provides a small perturbation into the solution and the simu- 
lation was run until time t = 1.5s. In Figure 2 the pressure perturbation and 
velocity profiles calculated based upon the second-order MF'W approach with 
400 numerical cells for the small disturbance were shown. Theses results have 
been also compared with the reference solution with 8000 numerical cells to- 
gether with the higher-order WENO scheme given in [9]. It is clear that the 
results are in good agreement with both reference solution and the higher- 
order WENO scheme, which confirms that even for a small-perturbation the 
effect of gravitational source term has been accurately treated within the 
flux-differences of the neighboring cells for the finite volume method. 


60 
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Figure 3: Numerical results obtained based on the MF'W approach for the small perturba- 
tion imposed for the polytropic steady-state solution in comparison with the higher-order 
WENO scheme and reference. solution (a) Pressure perturbation,(b) velocity 


Eventually the amount of A in the equation (@) was set to 0.1, which 
creates rather large amplitude disturbance into the solution. The simulation 
based on the second-order MFW approach with 200 computational cells and 
CFL=0.9 was then implemented until t=1.5s . Figure M shows plots for the 
pressure perturbations and velocity along with the results calculated using 
reference solution with 2000 cells and higher-order WENO scheme. As can be 
observed the large amplitude perturbation was also captured by the defined 
MFW approach and the results are in qualitative agreement with the higher- 
order WENO method. 


—*-MFW 
¢ Higher-Order-WENO 
——Reference Solution 


Pressure perturbation 


Figure 4: Numerical results obtained based on the MFW approach for the large perturba- 
tion imposed for the polytropic steady-state solution in comparison with the higher-order 
WENO scheme and reference. solution (a) Pressure perturbation, (b) velocity 
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Conclusions 


In this paper a modified flux-wave formula is presented for the solution one 
dimensional Euler equations including gravitational source term. The nu- 
merical solver defined herein is well-balanced and deals with gravitational 
source terms inside the differences between fluxes for the finite volume com- 
putational cells. In order to cope with difficulties with nonphysical results, 
a combination of exact and approximate Riemann solution was utilized. A 
modification for the Sod problem was first considered with the gravitational 
source term and the validations was made with the higher-order WENO 
scheme and rather identical results was achieved. Then, the well-balanced 
property of the scheme was tested and the ability of the method in model- 
ing the perturbation applied into the steady-state equilibrium in comparison 
with unbalanced scheme was demonstrated. For the problem with small and 
large amplitude of perturbations, the MFW approach gave the same results 
equal to higher-order WENO scheme and accurately captured disturbances 
imposed to the polytropic atmosphere. 
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